home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_4 / a4_1b.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.1 KB  |  127 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 4.1 (Evaluation of a Taylor Series).
  9. % Section    4.1,    Taylor Series and Calculation of Functions, Page 203
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program investigatges Taylor approximations.
  13.  
  14. %      Pn(x) = c(1) + c(2)x + c(2)x^2 + ... + c(n+1)x^n
  15.  
  16. % where the degree n of approximation is large  (n ~ 25).
  17.  
  18. % Coefficient lists for several functions have been
  19.  
  20. % stored in M-files named;   zcos  zsin  ztan  zexp
  21.  
  22. % zacos  zasin  zatan  zcosh   zsinh   zsqrt   zlog
  23.  
  24. % zsqrt4   zinv   zemx2d2   zcosde   zsinde   zlogq
  25.  
  26. pause    % Press any key continue.
  27.  
  28. clc;
  29. %       Approximations for  tan(x)
  30.  
  31. % Issue the command  ztan  to load the coefficients
  32.  
  33. % into the array  C.  The function name is loaded
  34.  
  35. % into the variable fun, the degree is loaded into  m.
  36.  
  37. % The endpoints of [a,b] are loaded into  a  and b.
  38.  
  39. % Load the Taylor coefficients.
  40.  
  41. [fun,dfun,ifun,x0,m,C,Ax] = ztan;
  42.  
  43. pause     % Press any key continue.
  44. clc;
  45. a = Ax(1,1);    % You can change the left  endpoint a.
  46. b = Ax(1,2);    % You can change the right endpoint b.
  47. c = Ax(1,3);
  48. d = Ax(1,4);
  49. n = m;          % You can change the degree n.
  50. C = flipud(C);
  51. C = C(1:n+1);
  52. C = flipud(C);
  53.  
  54. pause    % Press any key continue.
  55.  
  56. clc; clg;
  57. h = (b-a)/200;
  58. X = a:h:b;
  59. x = X; 
  60. Y = eval(fun);
  61. axis([a b c d]);...
  62. P = polyval(C,X);...
  63. plot(X,Y,'-g',X,P,'-r');...
  64. hold on;...
  65. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  66. xlabel('x');...
  67. ylabel('y');...
  68. Mx1 = ['Comparison of ',fun,' and P'];...
  69. Mx2 = [Mx1,num2str(n),'(x)'];...
  70. title(Mx2);...
  71. grid;...
  72. axis;...
  73. hold off;...
  74. shg; pause    % Press any key to continue.
  75.  
  76. Mx1 = 'The function is  f(x) = ';
  77. Mx2 = 'The interval is  ';
  78. Mx3 = 'Pn(x) = c(1)x^n + c(2)x^(n-1) + ... + c(n)x + c(n+1)';
  79. Mx4 = 'The degree is  n = ';
  80. Mx5 = ', and the coefficients list  C  is:';
  81. clc,format short e,echo off,diary output,...
  82. disp(''),disp([Mx1,fun]),...
  83. disp([Mx2,'[',num2str(a),' , ',num2str(b),']']),...
  84. disp(Mx3),disp([Mx4,num2str(n),Mx5]),...
  85. for i=1:5:n+1, disp(C([i:min(i+4,n+1)])'); end,...
  86. diary off, echo on
  87.  
  88. pause    % Press any key to graph f(x) - Pn(x).
  89.  
  90. clc; clg;
  91. a = -1;    % You can change the left  endpoint a.
  92. b =  1;    % You can change the right endpoint b.
  93. h = (b-a)/200;
  94. X = a:h:b;
  95. x = X; 
  96. Y = eval(fun);
  97. P = polyval(C,X);
  98. c = min(Y-P);
  99. d = max(Y-P);
  100. axis([a b c d]);...
  101. plot(X,Y-P,'-r');...
  102. hold on;...
  103. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  104. xlabel('x');...
  105. ylabel('y');...
  106. Mx1 = ['The error;  ',fun,' - P'];...
  107. Mx2 = [Mx1,num2str(n),'(x)'];...
  108. title(Mx2);...
  109. grid;...
  110. axis;...
  111. hold off;...
  112. shg; pause    % Press any key for a list of numerical computations.
  113.  
  114. clc; format long;
  115. a = -1.2;
  116. b =  1.2;
  117. X = a:0.1:b;
  118. x = X;
  119. Y = eval(fun);
  120. P = polyval(C,X);
  121. points = [X;Y;P;Y-P];
  122. Mx1=['Taylor polynomial approximation of f(x) = ',fun];
  123. Mx2='     x(k)               f(x(k))            Pn(x(k))           error';
  124. clc,echo off,diary output,...
  125. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  126. diary off,echo on
  127.